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Abstract 

In this work, we design an entropy stable, finite volume approximation for the shallow 
water magnetohydrodynamics (SWMHD) equations. The method is novel as we design an 
affordable analytical expression of the numerical interface flux function that exactly preserves 
the entropy, which is also the total energy for the SWMHD equations. To guarantee the 
discrete conservation of entropy requires a special treatment of a consistent source term for the 
SWMHD equations. With the goal of solving problems that may develop shocks, we determine 
a dissipation term to guarantee entropy stability for the numerical scheme. Numerical tests 
are performed to demonstrate the theoretical findings of entropy conservation and robustness. 


1 Introduction 

Systems of conservation laws that model a physical system or process often possess an additional 
conserved quantity, the entropy, which is conserved for smooth solutions but increases (or decreases 
according to the sign convention adopted) in the presence of shocks. A method is said to be entropy 
conservative if the local changes of entropy are the same as predicted by the entropy conservation 
law. The approximation is said to be entropy stable if it produces more entropy than an entropy 
conservative scheme. 

We introduce the important distinction between entropy conservation and stability because 
there is a problem with entropy conservative formulations. They may suffer breakdown if used 
without dissipation to capture shocks. Physically, entropy must be dissipated at a shock. However, 
an isentropic algorithm does not allow the capture of this physical process, which generates large 
amplitude oscillations around the shock [3]. A more dire issue is that entropy conservative formu¬ 
lations can not converge to the weak solution as there is no mechanism to admit the dissipation 
physically required at the shock. 

Tadmor [25, 26] was the first to introduce the idea of entropy conservation to design stable 
numerical approximations of nonlinear hyperbolic conservation laws. At a semi-discrete level the 
principle is that the discrete flux function satisfies discrete conservation laws of the conservative 
variables as well as the scalar conservation law for entropy in a discrete sense. Tadmor’s flux 
function performs well for smooth data but can be made stable for shock problems [27]. But, 
Tadmor’s flux function, which involves an integral in phase space, is numerically expensive. For 
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large scale simulations we want a more practical, computationally tractable, entropy conserving 
flux. 

Thus, in this paper, it is our goal to develop affordable, entropy stable methods for the shallow 
water magnetohydrodynamic (SWMHD) model. There is recent work for optimal, entropy stable 
approximations for the Euler equations from Ismail and Roe [15] which was extended to arbitrary 
order with a discontinuous Galerkin (DG) spectral element formulation by Garpenter et al. [3]. 
Also, there is recent work on entropy-stable, high-order, and well-balanced approximations for the 
shallow water equations [11, 12, 29]. 

The remainder of this paper is organized as follows: Sec. 2 provides an introduction to the 
SWMHD equations. In Sec. 3 we briefly describe the finite volume discretization used. Sec. 4 
defines the necessary variables and analytical tools to discuss the entropy in a mathematically 
rigorous way. We derive an entropy conserving numerical flux in Sec. 5.1 and discuss the stabilized 
flux in Sec. 5.2. Numerical results are presented in Sec. 6, where we verify theoretical predictions. 
Sec. 7 presents concluding remarks. Finally, in the Appendix we provide entropy stable fluxes that 
can be used for two dimensional computations. 


2 The Shallow Water Magnetohydrodynamic Model 

The shallow water magnetohydrodynamics (SWMHD) equations were first proposed by Gilman as a 
model for phenomena in the solar tachocline [13], which is the thin layer in a star between the outer 
turbulent convection zone, and the quiescent interior where heat transfer is predominantly radiative. 
The tachocline has many interesting physical applications such as the formation of sunspots [5]. 
Tachoclines can also serve as a useful model of solar dynamo action [7, 8, 13]. 

As was first shown by Gilman [13] with complete details presented by Rossmanith [24] the 
SWMHD equations are derived from the ideal MHD equations by integrating in the z—direction 
under the assumption of incompressibility, a two dimensional variation of the flow variables, and 
magnetohydrostatic equilibrium in the z—direction: 


dz 


p = constant, 


\B\ 


\B\ 


\B\\ 


= P9, 

= pgz, 


( 2 . 1 ) 


dz = pg 


Under these assumptions and the divergence-free constraint V • (hB) the SWMHD equations can be 
cast into the form of a system of conservation laws [5, 24]. The one dimensional SWMHD equations 
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comprise a hyperbolic system written in the conservative form 
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( 2 . 2 ) 


where the final equation is the one dimensional involution condition on the magnetic field B. The 
variables Vi, V 2 , and B in (2.2) are the horizontal components of the fluid velocity and magnetic 
field, h is the fluid layer depth, and g is the constant gravitational acceleration. 

Hyperbolic systems, like the SWMHD equations, are typically solved numerically with Godunov- 
type methods [20]. These methods require the solution of a Riemann problem at element interfaces. 
Unfortunately, the one-dimensional SWMHD equations are degenerate [5], just like the ideal MHD 
equations [22, 23], due to the involution, V • {hB) = 0. For example, if we consider the one 
dimensional problem we see from (2.2) that dt{hBi) = 0. Thus, there is no propagating Riemann 
wave associated with the conservative variable hBi. To alleviate such degeneracies we are interested 
in relaxing the involution V • (hB) — 0 such that a conventional Riemann solver based method may 
be used for the approximation of the SWMHD system. 

To accomplish this we weaken the divergence condition to require that V • (hB) « 0, i.e., the 
divergence condition is to be maintained up to finite precision error. Dellar [6] showed that for the 
relaxed assumption of the divergence condition we obtain a forced version of the SWMHD model 
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(2.3) 


where the source term is proportional to the divergence condition. Importantly, the source term 
in (2.3) alters the system to make hBi an advected scalar while maintaining the conservation of 
mass and momentum. The forced equations (2.3) can be derived from a Hamiltonian approach [6], 
so the system (2.3) conserves the total energy. The conservation of total energy is important when 
we develop the entropy conserving flux in Sec. 5.1. Also, the source term serves to restore Galilean 
invariance when V • (hB) yf 0 [6]. Finally, the source term in (2.3) is analogous to the Janhunen 
source term for the ideal MHD equations [16]. 

It is important to note that the source term in (2.3) preserves the conservation of the momentum 
and energy, while enforcing the involution condition on hB. Preserving the conservative properties 
of the fluid equations removes a significant drawback of the other common proposed Powell type 
source term for the SWMHD equations. The Powell source term renders the equations of momentum 
and energy for SWMHD non-conservative [6]. However, according to the Lax-Wendroff theorem [19] 
only conservative schemes can be expected to compute the correct jump conditions and propagation 
speed for a discontinuous solution. Difficulties of the non-conservative formulation obtaining the 
correct weak solution has been documented in the literature for the ideal MHD equations [28]. 
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3 Finite Volume Discretization 


The finite volume (FV) method is a discretization technique for partial differential equations espe¬ 
cially useful for the approximation of systems of conservations laws. The finite volume method is 
designed to approximate conservation laws in their original form, e.g., 


Utdx+ / f ■ hdS = 0. 


>dv 


(3.1) 


For instance, in one spatial dimension we break the interval into non-overlapping intervals (the 
“volumes”) 

r 1 

(3.2) 


X. 1jX. 1 
I— n 


2 ^2 

and the integral equation of a balance law, with a source term, becomes 

d 
dt 


=,+ 1/2 rXi+1/2 

udx + f* {xi+ 1 / 2 ) - f* {xi^i/ 2 ) = sdx 

- 1/2 ^^ i -1/2 


(3.3) 


At this point, we make approximations. A common approximation is to assume that the solution 
and the source term are constant within the volume. Then we determine, for example, what is 
analogous to a midpoint quadrature approximation of the solution integral 


i+1/2 


udx ; 
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Ui dx = UiAxi- 


(3.4) 


i-1/2 


i-1/2 


Note that the solution is typically discontinuous at the boundaries of the volumes. To resolve this, 
we introduce the idea of a “numerical flux, ” r{u\u^ ), often derived from the (approximate) 
solution of a Riemann problem. That is, f* is a function that takes the two states of the solution 
at an element interface and returns a single flux value. For consistency, we require that 


f*{u,u) = f, 


(3.5) 


that is, the numerical flux is equivalent to the physical flux if the states on each side of the interface 
are identical. A significant portion of this paper is devoted to the derivation of a numerical flux 
that conserves the discrete entropy of the system for the shallow water MHD equations. So we 
defer the discussion of the numerical flux to Sec. 5. 

We must also address how to discretize the source term s in (3.3). There is a significant amount of 
freedom in the source term discretization. But, previous work by Fjordholm et al. [9] demonstrated 
that designing entropy stable methods for the shallow water equations with discontinuous bottom 
topography required a special treatment of the source term. In the later derivations a consistent 
source term discretization necessary for entropy conservation will reveal itself. Thus, we also defer 
the discussion of the discrete treatment of the source term to Sec. 5.1. 


4 Entropy Analysis 

In this section we define the entropy variables and entropy Jacobian necessary to develop an entropy 
stable approximation for the SWMHD equations. We note that one can find a fully general and 
detailed description of entropy stability theory in, for example, [2, 10]. 
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For the SWMHD equations the total energy acts as an entropy function, [6] 


U (u) = i (^gh^ + hvi + hv2 + hB^ + hB^) , (4.1) 

where u is the vector of conserved variables. From the total entropy (4.1) we define the set of 
entropy variables 



The entropy variables (4.2) are equipped with the symmetric positive definite (s.p.d) Jacobian 
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(4.4) 


where = gh is the standard wave speed for the shallow water equations. The fact that the 
matrices H and are s.p.d. is a direct consequence that the entropy function (4.1) is a strongly 
convex, injective mapping [6]. For the derivations that follow we also require the flux Jacobian in 
the x—direction: 
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Because it will be of use in later derivations we also compute the entropy potential 

2 p = q ■ f - F = ^gh'^vi - hBi {viBi + V 2 B 2 ) , (4.6) 

where F is the entropy flux derived from the identity C/uA = F^, [9], 

F = gh^vi + ^ {hv\ + hviv^ + hviB"^ — hviBf) — hv 2 BiB 2 . (4.7) 
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4.1 Discrete Entropy Conservation for the ID SWMHD Equations 


We introduce the concept of discrete entropy conservation. Let’s assume that we have two adjacent 
states {L,R) with cell areas {Axl, Ax^). We discretize the SWMHD equations semi-discretely 
and examine the approximation at the t + | interface. We suppress the interface index unless it is 
necessary for clarity. 


Axl 


dUL 


Axr 


duR 

dt 


fL-f* + AxLSl, 
*+2 

f* - fn + AxRS i. 

*+2 


(4.8) 


We interpret the update (4.8) as a finite volume scheme where we have left and right cell-averaged 
values separated by a common flux interface. 

We premultiply the expressions (4.8) by the entropy variables to convert to entropy space. From 
the chain rule we know that Ut = q^Ut, hence a semi-discrete entropy update is 


Vi) ’ 


(4.9) 


If we denote the jump in a state as I'] = (•)ij — (-)L and the average of a state as= ((•)i^+(•)i)/2, 
then we find the total element update to the entropy is 

^^{AxrUl + AxrUr) = - |q • /] 4- {qf' f* + 2 {AxqJ s^^i. (4.10) 

We want the discrete entropy update to satisfy the discrete entropy conservation law. To achieve 
this, we require 

2 fAxqJ^ S i -{q- fl + f* = - [T’l ■ (4.11) 

We combine the known entropy potential ip in (4.6) and the linearity of the jump operator to rewrite 
the entropy conservation condition (4.11) as 

f* {viBi + V 2 B 2 )j - 2 {Axq}'^ s^^. (4.12) 

We denote the constraint (4.12) as the discrete entropy conserving condition. However, this is a 
single condition on the vector /*, so there are many potential solutions for the entropy conserving 
flux. However, we have the additional requirement that the numerical flux must be consistent (3.5). 
We develop an affordable entropy conserving flux function in Sec. 5. We also defer the discretization 
of the source term Si^i /2 to to the next section because special care must be taken to guarantee 
that (4.12) is satisfied. 


5 Derivation of an Entropy Stable Numerical Flux 

With the necessary entropy variable and Jacobian definitions as well as the formulation of the dis¬ 
crete entropy conserving condition (4.12) we are ready to derive an affordable entropy conserving 
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numerical flux in Sec. 5.1. As we previously noted, entropy conserving methods may suffer break¬ 
down in the presence of shocks [3]. Thus, in Sec. 5.2, we will design an entropy stable Riemann 
solver that uses the entropy conserving flux as a base and incorporates a dissipation term required 
for stability. 


5.1 Entropy Conserving Numerical Flux for ID SWMHD 

In this section we define an affordable, entropy conserving numerical flux for the one dimensional 
SWMHD equations. 


Theorem 1. (Entropy Conserving Numerical Flux) If we discretize the source term in the finite 
volume method to contribute to each element as 
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then we can determine a discrete, entropy conservative flux to he 
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Proof. To derive an affordable entropy conservative flux for the one dimensional SWMHD equations 
we first expand the discrete entropy conserving condition (4.12) componentwise to And 
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To determine the unknown flux components of f* we use (5.3) to obtain a system of equations for 
each of the linear jump terms. For this expansion step we make repeated use of the linear jump 
operator properties: 

[a61 = m + m H , = 2 laj [al. (5.4) 

We expand the terms of (5.3) to linear jump terms of the primitive variables h, vi, V 2 , Bi, and H2- 
This is straightforward in all but the last two terms on the right hand side 

9 w /r + (/2 - /r) bii + if; - s^^2:& /d + u: - /d ii?ii + u; - ib^j /d m 

= {b^}} + 9 M - IhBiiviBi + V 2 B 2 )j - 2 , i • 

(5.5) 
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Care must be taken when we expand the the second term on the right hand side of the entropy 
condition (5.5). If we naively split the jump terms all the way to individual linear jump terms it will 
introduce averages of the magnetic field variables B on the jump in fluid height h. However, the 
continuity equation in the SWMHD model contains no influence from the magnetic held. Therefore, 
it would be impossible for the resulting numerical flux to be consistent with the physical flux. 

To remove this inconsistency we split the magnetic held jumps on the right hand side of (5.5) 
in such a way that cancellation with the source term is possible. Thus, rather than fully expanding 
each jump term down to the primitive variable level, it is useful to leave the product hBi term as 
a jump. This choice is motivated by the structure of the source term because 
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We expand the second term on the right hand side of (5.5) as follows: 

lhBi{viBi + U2H2)] = + I/1U2.B1H2], 

= ihB^J IvM + ivM IhB,} + ihB^J [U2B2I + IV 2 B 2 } IhB,} , 

= (IvM + iv^B,} ) lhB,j + ihB,} [uiHil + ihB,} [U2H2I , (5.7) 

= ( iviBil + iv2B2l ) I/rHil + IhB.l ( ^ {B,} + [uj ) 

+ ihB,}{iv2}lB2} + lB2}lv2}). 


We substitute (5.7) into the expanded entropy constraint (5.5) to obtain 

g w /r + (/2 - fD bii + ifs - /r) + (/; - /r) m + in - ib,} /*) m 

= Ig {{h^ +g m 1/^1 - ihB^} (m + ib,} [uii + su 2 :& [B2I + ib,} [U2I) 


- ( iviBil + iv 2 B 2 } ) IhB.l - 2 


i-\-7 


(5.8) 

From the structure of the expanded entropy conserving constraint (5.8) we are prepared to select 
a discretization of the source term at the interface. The goal is to cancel the jump in hBi term. 
We choose that the source term discretization on cell i will contribute to the interfaces i + ^ and 
i — h- In exact arithmetic the dot product of the entropy variables with the source term is 


- (^1^1 + ^2-82). 


(5.9) 


From (5.9) we see that it will be possible to create a consistent discretization to cancel the 
term that remains in (5.8). To do so we require the source term discretization 
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The source term discretization (5.10) is consistent to the source term in (2.3). For example, on a 
regular grid where Axl = Axr we have for Si^i /2 
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where consistency is easily checked. We note that there may be computational issues if either Bi 
or B 2 is zero. However, if i ?2 = 0 (for example), then the contribution to the entropy conservation 
constraint is also zero. Now it is possible to use the source term discretization (5.10) to see the 
explicit cancellation of the problematic term at the i + ^ interface: 


-iiviBiJ + iv2B2} ) lhB,l - 2 SAxqf 


s^^=-{ivM + iv2B2l)lhB^l 


-2lAxq} 


lhB,l 


V 


0 1\ 
0 
0 

I f^>2B2}l I I 

LgAa:B2}l J / 


= -{lv,B,l + iv2B2l)lhB,l 
+ {ivM + iv2B2})lhB,j, 

= 0 . 

(5.12) 

With the remaining |hi?i] term is removed from (5.8) we are left with the constraint on the 
numerical flux that only depends on the linear jump in each primitive variable 


9 w /r + (/2* - /D bii + ifs - /r) i«2i + in - /r) i^^ii + (/s - id m 

= I9 {m} i^'ii+5 Mm w - ihB^} (m + iB^} [uii + su 2 :& m + m 1^21). 

(5.13) 

We collect the like terms and determine a set of equations for the unknown components of the 
entropy conserving numerical flux 
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The equations (5.14) are easily solved to determine the form of the entropy conserving flux to be 
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Interestingly, to guarantee the conservation of entropy we had to relax the divergence-free condition 
to dx{hBi) « 0. This assumption manifests itself in the fourth component of which is nonzero. 
Essentially the fourth term in is the commutator of multiplication and averaging of hBi. ■ 


Remark 1. If the flow occurs without the presence of a magnetic held then the SWMHD model be¬ 
comes the traditional shallow water (SW) equations. The structure of /*>'='= (5.15) can be separated 
into a shallow water component and magnetic held component, i.e.. 
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and becomes the appropriate entropy conserving flux for the SW equations described in [9, 12]. 

Remark 2. We have considered a computational problem with a flat bottom topography. So it is 
a natural question for a shallow water type model if the approximation remains well-balanced, i.e., 
recovers a constant fluid height solution in the presence of a non-constant bottom topography [29]. 
Previous work by Fjordholm et al. [9] provides a consistent, well-balanced, and entropy conserving 
treatment of a non-constant (and even discontinuous) bottom topography b. For example in the 
a;—direction 



In the presence of a non-constant bottom topography one would add the discretized source term 
(5.18) to the X component of the momentum equation of the SWMHD model. Note that the 
treatment of the bottom topography source term by Fjordholm et al. in (5.18) is analogous to our 
treatment of the Janhunen source term for SWMHD (5.10). 
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5.2 Dissipation Terms for an Entropy Stable Flnx 

To create an entropy stable numerical flux function we use the entropy conserving flux (5.15) as a 
base and subtract some form of numerical dissipation, e.g, 

r = r’“-^DH, (5.19) 

where D is a dissipation matrix. 

Our goal to guarantee the entropy stability of the approximation is to reformulate the dissipation 
term (5.19) such that it is guaranteed to be positive. To do so we will consider the eigenstructure 
of the flux Jacobian A altered by the Powell source term 
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It is important to note that we use the Janhunen source term (2.3) to derive an entropy conservative 
numerical flux function in Sec. 5.1. However, to design an entropy stable approximation we require 
that the eigendecomposition of the flux Jacobian matrix can be related to the entropy Jacobian 
(4.4). This particular scaling, first examined by Merriam [21] and explored more thoroughly by 
Barth [2], requires that the PDE system is symmetrizable. Previous analysis of the SWMHD system 
by Dellar [6] and ideal MHD equations [2, 14] have demonstrated that the Powell source term is 
necessary to restore a symmetric MHD system. We reiterate that the altered flux Jacobian is used 
only to derive the dissipation term. Just as Lax-Friedrichs differs from Roe in the construction 
of a dissipation term, we use the Powell source term only to build our dissipation term. Thus, 
we do not introduce any inconsistency in the previous entropy conserving flux derivations. We 
demonstrate with numerical tests in Sec. 6.4 that the newly developed entropy stable flux functions 
are competitive with the well known Roe flux. 

We denote the flux Jacobian modified with the Powell source term (5.20) as 
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equipped with a full set of eigenvalues 


Ai=vi-Cg, X2 = vi-Bi, A3 = ui, A4 = ui + i?i, A5 = ui + Cg, (5.22) 
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where = gh + is the magnetogravity wave speed. 

Now that we have a symmetrizable matrix A we utilize a previous result from Barth [2] which 
provides a systematic approach to restructure a general eigenvalue problem to a symmetric eigen¬ 
value problem. To do so we rescale the right eigenvectors of an eigendecomposition with respect to 
a right symmetrizer matrix in the following way: 

Lemma 1. (Eigenvector Scaling) Let A e he an arbitrary diagonalizahle matrix and S the 

set of all right symmetrizers: 

5= Biss.p.d, AB={ABf]. (5.24) 

Further, let R G denote the right eigenvector matrix which diagonalizes A, i.e., A ~ RAR~^, 

with r distinct eigenvalues. Then for each B G S there exists a symmetric block diagonal matrix T 
that block scales columns of R, R = RT, such that 

B = RR^, A = RAR~\ (5.25) 

which implies 

AB=RAIt. (5.26) 

Proof. The proof of the eigenvector scaling lemma is given in [2]. ■ 

Theorem 2. (Entropy Stable 1 (ESI)) If we apply the diagonal scaling matrix 
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to the matrix of right eigenvectors R (5.23), then we obtain the Merriam identity [21] (Eq. 
P9- 77 ) 
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that relates the right eigenvectors of A to the entropy Jacobian matrix (4.4). For convenience, we 
introduce the diagonal scaling matrix S = in (5.28). We then have the guaranteed entropy stable 
flux interface contribution 
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Proof. We define the dissipation term in the numerical flux (5.19) to be 
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where the eigendecomposition of A = RAL is given by (5.22) and (5.23). We define entropy 
stability to mean the approximation guarantees that the entropy within the system is a decreasing 
function, satisfying the following inequality 


dU 

~di 


dF 

dx 


q^s < 0. 
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(5.31) 




From the previously computed discrete entropy update (4.10) the total entropy within an element 
(now including the dissipative term (5.30)) is 

^{AxlUl + AxrUr) = - |q • /] + f* + 2-8;AxqS- , 

= -Iq- fl + yf /*’“ - ^ R| A|L |m] + 2 lAxqS- i, 

+ AxrUr) + m=-^ Iqf R|A|L M , 

(5.32) 

from the design of the entropy conserving flux To ensure entropy stability, we must guarantee 

that the RHS term in (5.32) is non-positive. Unfortunately, it was shown by Barth [2] that the 
term 

-i[qfR|A|LM, (5.33) 

may become positive in the presence of very strong shocks. However, we know from entropy sym- 
metrization theory, e.g [2, 21], that the entropy Jacobian H, given by (4.4), is a right symmetrizer 
for the flux Jacobian that incorporates the Powell source term A. Therefore, with the proper scaling 
matrix T we acquire the Merriam identity 

H = RR^ = (rt) (rt)^ = RSR^. (5.34) 


The rescaling of the right eigenvectors of A to satisfy the Merriam identity (5.28) is sufficient to 
guarantee the negativity of (5.33). We see from (5.33) and (5.34) 

^ |q]^R|A|LMq |q] , 

= -i[qf R|A|LH[ql, 

f ^ (5-35) 

= --[qf R|A|L(feR^) [ql, 

= -i[qf R|A|SR^ [ql, 

^ ^ -1 

where we used that L = R .So with the appropriate diagonal scaling matrix S we have shown 
that (5.35) is guaranteed negative because the product is a quadratic form scaled by a negative. 
We use the right eigenvectors from (5.23), the constraint (5.28), and straightforward algebraic 
manipulation to determine the diagonal scaling matrix 


T = diag 


c c 
Cgy/2g ’ y2q ’ 


Bi c c \ 
Cgy/g ’ ’ Cgy/^j ’ 


(5.36) 


where = gh and = gh+Bf are the wave speed and magnetogravity wave speed respectively. I 

Remark 3. There are other possible, negativity guaranteeing (but more dissipative) choices for the 
dissipation term (5.19). For example, if we make the simple choice of dissipation matrix to be 




(5.37) 


13 




where Xmax is the largest eigenvalue of the system from (5.22) and I is the identity matrix, then 
we obtain a local Lax-Friedrichs type interface stabilization 


= r’“-^|A™a.|H[gl, (5.38) 

= r’--^R|A^„,|SR^[ql, 

where, again, we use the Merriam identity (5.28) for the entropy Jacobian H. 


6 Numerical Results 

In this section, we numerically verify the theoretical findings for the entropy conserving and entropy 
stable approximations for the SWMHD equations. To integrate the semi-discrete formulation in 
time we use a low storage five-stage, fourth-order accurate Runge-Kutta time integrator of Car¬ 
penter and Kennedy [4]. First, in Sec. 6.1, we consider a test problem with a known analytical 
solution to demonstrate the accuracy of the entropy conserving method as well as the optimal and 
Lax-Friedrichs stabilized formulations. Next, Sec. 6.2 demonstrates the entropy conservation of the 
approximation for a strong Riemann problem. In Sec. 6.3 we demonstrate the computed solution of 
a strong Riemann problem for the entropy conserving method. This solution will exhibit significant 
oscillations in shocked regions. Finally, in Sec. 6.4 we compare the two entropy stable approxi¬ 
mations (5.29) and (5.38) against a high-resolution approximation comparable to that presented in 
the literature [5, 24]. 


6.1 Convergence 

For the convergence test, we switch to the manufactured solution technique and generate a smooth 
and periodic solution 


h{x, t) 


2 -1- sin (27r {x — t)) 

hvi{x,t) 


2 -1- sin (27r (x — t)) 

hv2{x,t) 

= 

2 -1- sin (27r (x — t)) 

hBi{x, t) 


1 

hB2\x, t) 


4 -1- 2 sin (27r (x — t)) 


with an additional analytic source term on the right hand side 


si{x,t) 


0 

/ \ 

S2{x,t) 


hx{x,t) [gh{x,t) + 

S3{x,t) 

= 

0 

S4(x,<) 


0 

S 5 {x,t)_ 


0 


on the domain = [—1,1] with periodic boundary conditions and final time T = 2. We select the 
time step for the RK method small enough such that the error in the approximation is dominated 
by the error in the spatial discretizations. For all computations, a regular grid is chosen according 
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to the number of grid cells listed in the tables below. Except for the entropy conserving scheme, 
where we also test a stretched mesh with a fixed ratio of 


AXr, 


AXr 


= 4. 


(6.3) 


First, we test the entropy conserving scheme on a regular grid. The L 2 -errors for all conserved 
quantities are shown in Tbl. 1. It is interesting to note, that for the specific regular grid used in 
the numerical experiment, the entropy conserving scheme achieves second order accuracy. We note 
that the average accuracy is hovering around 1.7, however looking closely at the finest grid results, 
it is clear that the experimental order of convergence is close to 2. However, the higher order 
convergence for the finite volume scheme is an effect of approximating the solution on a regular 
grid. The scheme drops to first order accuracy when an irregular grid is chosen as shown in Tbl. 2, 
where we stretch the grid with a constant factor of (6.3). 


# elements 

L 2 error h 

L 2 error hvi 

L 2 error hv 2 

L 2 error hBi 

L 2 error hB 2 

50 

2.05E-2 

2.12E-2 

2.91E-2 

3.77E-2 

1.62E-1 

100 

8.78E-3 

9.14E-3 

9.45E-3 

1.90E-2 

3.08E-2 

200 

2.45E-3 

2.50E-3 

2.59E-3 

3.12E-3 

7.22E-3 

400 

6.17E-4 

6.26E-4 

6.37E-4 

7.85E-4 

1.79E-3 

avg EOC 

1.69 

1.69 

1.84 

1.86 

2.17 


Table 1: Experimental order of convergence EOC for entropy conserving scheme on a regular grid 


# elements 

L 2 error h 

L 2 error hvi 

L 2 error hv 2 

L 2 error hBi 

L 2 error hB 2 

50 

1.44E-1 

2.50E-1 

5.90E-1 

1.64E-1 

1.05E-0 

100 

5.57E-2 

9.66E-2 

2.11E-1 

5.57E-2 

2.62E-1 

200 

2.73E-2 

4.57E-2 

1.02E-2 

2.47E-2 

1.22E-1 

400 

1.34E-2 

2.26E-2 

5.45E-2 

1.20E-2 

6.22E-2 

avg EOC 

1.14 

1.16 

1.15 

1.26 

1.36 


Table 2: Experimental order of convergence EOC for entropy conserving scheme on an irregular 
grid with a stretching factor of (6.3). 

Next, we demonstrate the convergence of the two entropy stable finite volume schemes. The 
convergence results for the ESI method are shown for the regular grid test in Tbl. 3 and the irregular 
grid test in Tbl. 4, where we see the average experimental convergence rate for either the regular or 
irregular grids are both first order accurate. The convergence behavior for the ES2 scheme on both 
grids fluctuates. Thus, we decided to run more tests on even finer grids to demonstrate that, again, 
the experimental order of convergence for regular or irregular grids is in the range of 0.7 — 1.36. 

6.2 Mass, Momentum, and Entropy Conservation 

By design we know the finite volume scheme for the SWMHD equations with the Janhunen source 
term (2.3) will conserve mass and momentum. From the derivations in 5.1 we know that the scheme 
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# elements 

L 2 error h 

L 2 error hvi 

L 2 error hv 2 

L 2 error hBi 

L 2 error hB 2 

50 

6.24E-2 

1.87E-1 

6.44E-2 

6.84E-3 

1.16E-1 

100 

2.91E-2 

8.49E-2 

2.94E-2 

3.20E-3 

5.33E-2 

200 

1.25E-2 

3.38E-2 

1.25E-2 

1.33E-3 

2.28E-2 

400 

5.79E-3 

1.27E-2 

5.78E-3 

4.72E-4 

l.lOE-2 

avg EOC 

1.14 

1.29 

1.16 

1.29 

1.13 


Table 3: Experimental order of convergence EOC for ESI scheme (5.29) on a regular grid. 


# elements 

L 2 error h 

L 2 error hvi 

L 2 error hv 2 

L 2 error hBi 

L 2 error hB 2 

50 

8.29E-2 

2.72E-1 

8.81E-2 

9.06E-3 

1.66E-1 

100 

4.25E-2 

1.40E-1 

4.40E-2 

4.20E-3 

8.33E-2 

200 

2.04E-2 

6.66E-2 

2.07E-2 

1.82E-3 

3.96E-2 

400 

1.02E-2 

3.147E-2 

1.02E-2 

6.66E-4 

1.99E-2 

avg EOC 

1.01 

1.04 

1.04 

1.26 

1.02 


Table 4: Experimental order of convergence EOC for ESI scheme (5.29) on an irregular grid with 
a stretching factor of (6.3). 


# elements 

L 2 error h 

L 2 error hvi 

L 2 error hv 2 

L 2 error hBi 

L 2 error hB 2 

50 

2.29E-2 

2.19E-1 

2.48E-2 

2.91E-3 

4.80E-2 

100 

1.02E-2 

9.92E-2 

1.14E-2 

1.41E-3 

2.13E-2 

200 

9.14E-3 

4.03E-2 

9.37E-3 

6.74E-4 

1.85E-2 

400 

7.06E-3 

1.62E-2 

7.08E-3 

2.92E-4 

1.42E-2 

800 

4.44E-3 

7.06E-3 

4.44E-3 

1.08E-4 

8.89E-3 

1600 

2.49E-3 

3.32E-3 

2.49E-3 

3.48E-5 

4.99E-3 

3200 

1.32E-3 

1.62E-3 

1.32E-3 

l.OlE-5 

2.64E-3 

avg EOC 

1.24 

1.18 

0.70 

1.36 

0.70 


Table 5: Experimental order of convergence EOC for ES2 scheme (5.38) on a regular grid. 


also exactly preserves the entropy on a general grid, if we use the newly designed flux (5.15). 
We measure the change in any of the conservative variables with 

Ae(r) := \eint{t = 0) - e„t(T)|, (6.4) 

where, for example, Cmt is the entropy function 

^ ~ \ ^^2) (6-5) 

integrated with a Gauss-Lobatto quadrature rule over the whole domain. 

To demonstrate the conservative properties we consider the one dimensional strong Riemann 
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# elements 

L 2 error h 

L 2 error hvi 

L 2 error hv 2 

L 2 error hBi 

L 2 error hB 2 

50 

6.65E-2 

3.13E-1 

6.74E-2 

3.66E-3 

1.38E-1 

100 

3.39E-2 

1.58E-1 

3.43E-2 

1.81E-3 

6.98E-2 

200 

1.97E-2 

7.29E-2 

1.97E-2 

8.82E-4 

3.96E-2 

400 

1.26E-2 

3.37E-2 

1.26E-2 

3.94E-4 

2.51E-2 

800 

7.83E-3 

1.64E-2 

7.80E-3 

1.51E-4 

1.56E-2 

1600 

4.57E-3 

8.22E-3 

4.57E-3 

4.97E-5 

9.12E-3 

3200 

2.54E-3 

4.15E-3 

2.54E-3 

1.48E-5 

5.07E-3 

avg EOC 

0.78 

1.04 

0.79 

1.33 

0.79 


Table 6: Experimental order of convergence EOC for ES2 scheme on an irregular grid with a 
stretching factor of (6.3). 


problem taken from [5] with the initial condition 


h 

Vl 

V2 

Bi 

B2 


[1,0,0, if a;<0, 

[2,0,0,0.5,1]^, if coO, 


( 6 . 6 ) 


on the domain = [—1,1] and periodic boundary conditions. We compute the error in the change 
of each conservative variable for three values of the Courant-Friedrichs-Lewy (CEL) number: 1.0, 
0.1, and 0.01. For each simulation the entropy conserving finite volume scheme used 100 regular 
grid cells for the results in Tbl. 7 and 100 irregular grid cells in Tbl. 8. The final time of each 
simulation was T = 0.4. 

We see that for each value of the CEL number we obtain errors in the mass and momentum 
on the order of machine precision. Also, as expected, the addition of the Janhunen source term to 
enforce the divergence condition causes us to lose conservation of the quantities hB. Interestingly, 
in the error for the change in entropy we see the dissipative influence of the time integrator. The 
dissipation introduced by the temporal discretization can be significantly reduced if we shrink the 
time step. We see that the error of the entropy can be lowered to single or double machine precision 
by decreasing the CEL number, and hence the time step. 

It is also possible to see the temporal accuracy in Tbl. 7 and Tbl. 8. If we shrink the time step 
by a factor ten we see that the error in the entropy shrinks by a little over a factor of 10"^, as we 
expect for a fourth order time integrator. 


6.3 Entropy Conserving Riemann Problem 

We now apply the entropy conserving (EC) scheme to the strong Riemann problem (6.6), except 
we choose inflow/outflow type boundary conditions. The computation is performed on 100 regular 
grid cells with CEL number 0.1 up to a final time of T = 0.4. To create a reference solution we use 
the ESI scheme on 5000 grid cells. We present the solution computed with the EC finite volume 
method against the reference solution in Fig. 1. 

The results for the computed fluid height h, velocities vi, V 2 , and magnetic field Bi, B 2 respec¬ 
tively in Fig. 1 show that the EC scheme computes the rarefactions and the shocks in the solution 
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CFL 

h 

hvi 

hV2 

hBi 

hB2 

U 

1.0 

-2.22E- 16 

3.53F;- 14 

-5.06E- 14 

-2.82E - 04 

6.61F;- 04 

-7.88E - 05 

0.1 

4.44F; - 16 

8.28E- 14 

-6.99E- 14 

-2.82E - 04 

6.64F; - 04 

-1.54E - 09 

0.01 

-9.99F;- 16 

9.66F;- 14 

-7.48F; - 14 

-2.82E - 04 

6.64F; - 04 

-9.20E - 14 


Table 7: Conservation errors (integrated over the whole domain) of the entropy conserving approx¬ 
imation of the Riemann problem (6.6) for different CFL numbers, final time T = 0.4, and 100 
regular grid cells. 


CFL 

h 

h vi 

hV2 

hBi 

hB2 

U 

1.0 

0.0 

0.0 

4.51E- 17 

-9.67E - 04 

2A7E - 03 

-2.05F;- 05 

0.1 

-3.33E - 16 

-1.04E- 17 

-2.57E - 14 

-9.66F;- 04 

2.47F; - 03 

-3.74F; - 10 

0.01 

-3.44E - 15 

3.40E - 16 

2.64E- 16 

-9.66F;- 04 

2.47F; - 03 

-2.HE - 14 


Table 8: Conservation errors (integrated over the whole domain) for different CFL numbers with 
end time t = 0.4 on an irregular grid with a stretching factor of (6.3). 


quite accurately, but at the expense of large post-shock oscillations. These oscillations are expected 
as energy must be dissipated across the shock but the EC scheme is basically dissipation free except 
for the influence of the time integrator (as we previously noted in Sec. 6.2). 

6.4 Entropy Stable Riemann Problem 

Next, we use the ESI and ES2 schemes to compute the solution of (6.6) with inflow/outflow type 
boundary conditions. Again, the computation is performed on 100 regular grid cells with CFL 
number 0.1 up to a final time of T = 0.4 with a reference solution created using a high-resolution 
run of the ESI scheme on 5000 grid cells. We note that the same reference solution could be obtained 
from a high-resolution computation of the ES2 scheme. We present the solution computed with the 
entropy stable methods against the reference solution in Fig. 2. 

Our computed results for the computed fluid height h, velocities ui, V 2 , and magnetic field Bi, 
i ?2 shown in Fig. 2 compare well to the literature [5, 24]. We see that, as was discussed in Sec. 5.2, 
the Lax-Friedrichs type stabilization introduces more dissipation into the approximation. 

As a final test we compare the ESI scheme to a standard Roe scheme. The development of 
Roe type schemes for the SWMHD equations has been considered by several authors and complete 
details can be found in [5, 17, 24]. The comparative computation used 100 regular grid cells with 
a CFL number of 0.1 up integrated to T = 0.4. A reference solution was created from a high- 
resolution run of the ESI scheme on 5000 grid cells. The computed results of the schemes are 
presented in Fig. 3. We see that for the same number of degrees of freedom the ESI scheme is less 
dissipative than the Roe scheme. This is particularly evident in the approximation of quantities in 
the tangential direction, V 2 and i? 2 - We also note that the Roe scheme is more dissipative than the 
ES2 Lax-Friedrichs type entropy stable approximation which can be seen from a comparison of the 
results in Fig. 2 and Fig. 3. 
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7 Concluding Remarks 


In this work we present a novel, affordable, entropy stable ffux for the shallow water MHD equations. 
The derivations required us to relax the involution condition such that dx{hBi) « 0 within finite 
precision error. Under this less restrictive assumption we were able to derive an entropy conserving 
flux, denoted Special care had to be taken for the discretization of the source term in order 

to guarantee discrete entropy conservation. Because entropy conserving approximations can suffer 
breakdown at shocks we extended our analysis and derived two stabilizing dissipation terms that 
we apply to the entropy conserving flux. We finally used a variety of numerical test examples to 
demonstrate and underline the theoretical findings. 

The derivation of the entropy conserving and entropy stable numerical fluxes in this paper fo¬ 
cused on the one-dimensional SWMHD equations. The restriction to one spatial dimension was 
because the analysis proved to be quite involved. However, the discussion was general, so the 
derivations in this paper readily extend to provably entropy conserving and entropy stable approx¬ 
imations of the SWMHD equations on multi-dimensional Cartesian grids. In the Appendix we 
provide details on the multidimensional flux functions. 

It is well known that the issue of divergence-free constraint and the accuracy of an approximation 
is more serious in two dimensions. As was discussed in Sec. 2 the Janhunen source term acts, 
analogously, to a hyperbolic divergence cleaning method where the error in the involution term is 
radiated out of the computational domain. To verify the utility of our method we consider a two 
dimensional example similar to the “rotor problem” of Balsara and Spicer [1] and considered by 
Toth [28] for the ideal MHD equations. The computational domain is H = [— 1,1] X [-1, 1] with 
initial data 

\hA 


Vi 


V2 

Bi 

B2 


[10,-j/, a:, 0.1,0]^, if ||a;||2<0.1, 
[1,0,0,1,0]^, if ||a;||2>0.1, 


(7.1) 


integrated up to the final time T = 0.2. The computational results of the two dimensional ESI 
scheme with 200 regular grid cells in each direction are presented in Fig. 4. We note that the 
computed fluid height h and velocity components vi and V 2 compare well with previous results 
found in the literature [18]. As we expect, because we use a different divergence cleaning procedure, 
our results for the Bi and B 2 magnetic fields differ slightly from those computed in [18]. 

In this paper we have demonstrated that the entropy stable approximations for the SWMHD 
equations are competitive with existing solvers. We also showed that the methods readily extend 
to multiple dimensions because the source term we need for entropy conservation also acts as a 
divergence cleaning technique. 


A Flux Functions in Two Spatial Dimensions 

We note that, for the discussion of the two dimension fluxes we suppress the i index on the multi¬ 
dimensional approximation for the purpose of clarity. 
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A.l Entropy Conservative Flux 

To develop the entropy conservative flux in the y—direction we first note that the entropy potential 

-tpy = V g - G = - hB 2 {viBi + V 2 B 2 ), (A.l) 

where we have the entropy flux in the y—direction 

G = gh^V2 + ^ {hv‘lv2 + hv^ + hv2Bf — hv2B2) — hviBiB2. (A.2) 


Note that the discrete entropy conservation condition (4.11) has the same structure in each Carte¬ 
sian direction. Lastly, the source term contributes symmetrically to each direction. With a proof 
analogous of that for in Sec. 5 we present the entropy conserving numerical flux for the 

y—direction. 


Corollary 1. (Entropy Conserving Numerical Flux: y—direction) If we discretize the source term 
in the finite volume method to contribute to each element as 


s 


j ~ 


1 

2 



/ 

0 


0 

\ 


0 


0 


2 

0 

+ lhB2l^ 1 

J 2 

0 

f^iSiS 


V 

fiAxBiS 

fv2B2} 

L«AxB2»J 


iv2B2'^ 

t-h/ 


then we can determine a discrete, entropy conservative flux of the form 


g 


*,ec 


mM 

mivi}iv2}~ihB2jiB,J 

miv2r+h{h^]-ihB2}iB2} 

miv2}iB^J-ihB2jM 

miv2}iB2}-ihB2}{{v2} 


(A.4) 


A.2 Entropy Stable Fluxes 

Just as was done in Sec. 5.2 we can create 2D entropy stable flux functions. We, again, use the 
flux Jacobian altered by the Powell source term to create the dissipation term: 


0 

0 

1 

0 

0 

— VxV2 + B 1 B 2 

V2 

Vl 

-B 2 

0 

Cq 

-f 

1 

0 

2v2 

0 

-B- 

VIB 2 — V 2 BI 

-B 2 

Bi 

V2 

0 

0 

0 

0 

0 

V2 


equipped with a full set of eigenvalues 

Al = W2 - Cg, \2=V2- B 2 , A3 = V2, A4 = W2 + B 2 , A5 = U 2 + Cg, 


(A.5) 


(A.6) 
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and right eigenvectors 



■ 1 

0 

1 

0 

1 


Vi 

1 

Vl 

1 

Vl 

Rj, — 

V2 - Cg 

0 

V2 

0 

V2 + C, 


Bi 

1 

Bi 

-1 

Bi 


0 

0 

4 

B2 

0 

0 


(A.7) 


where = gh + is the magnetogravity wave speed. 

Corollary 2. (Entropy Stable 1 (ESI): y—direction) If we apply the diagonal scaling matrix 


Ty = diag 


c c 

Cgy/2g ’ ./Ig 


B 2 c c \ 


(A.8) 


to the matrix of right eigenvectors Ry (A.7), then we obtain the Merriam identity [21] (Eg. 7.3.1 
pg- 77) 

H= RyRy = (^RyTy^ ^Ry Ty^ = RySyRy , (-^- 9 ) 

that relates the right eigenvectors of B to the entropy Jacobian matrix (4.4). For convenience, we 
introduce the diagonal scaling matrix Sy = Ty in (A.9). We then have the guaranteed entropy stable 
flux interface contribution 

g*,ESl ^ g*,ec _ 1 ^ g*,ec _ |qj , (A.IQ) 

Remark 4. (Entropy Stable 2 (ES2): y—direction) If we choose the dissipation matrix to be 


D= |A, 


(A.ll) 


where \max is the largest eigenvalue of the system from B and I is the identity matrix, then we 
obtain a local Lax-Friedrichs type interface stabilization 


g*,ES 2 _ 




|HM. 


(A.12) 
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Figure 1: In blue we present the entropy conserving approximations of the quantities h, vi, V 2 , -Bi, 
and B 2 at T = 0.4. Solid black represents the reference solution. 











































Figure 2: The entropy stable approximations of the quantities h, vi, V 2 , Bi, and i ?2 at T = 0.4. 
Solid black is the reference solution, dashed blue is the ES2 scheme, and red with knots is the ESI 
scheme. 







































Figure 3: The entropy stable results compared against a Roe scheme for the quantities h, vi, V 2 , 
Bi, and i ?2 at T = 0.4. Solid black is the reference solution, dashed blue is the Roe scheme, and 
red with knots is the ESI scheme. 
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Figure 4: The computed solution to the “rotor problem” (7.1) using the ESI scheme for the quan¬ 
tities h, vij U 2 , Bi, and .S 2 at T = 0.2. Dashed lines are negative contour lines. 


























